Universal SSE algorithm for Heisenberg model and Bose Hubbard model with 

interaction 
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We propose universal SSE method for simulation of Heisenberg model with arbitrary spin and 
Bose Hubbard model with interaction. We report on the first calculations of soft-core bosons with 
interaction by the SSE method. Moreover we develop a simple procedure for increase efficiency of 
the algorithm. From calculation of integrated autocorrelation times we conclude that the method 
is efficient for both models and essentially eliminates the critical slowing down problem. 
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I. INTRODUCTION 

Recently, significant progress in quantum Monte Carlo 
methods has been observed. During the last two decades, 
advanced quantum Monte Carlo algorithms have been de- 
veloped. First quantum Monte Carlo methods, so-called 
world line algorithms, were based on Suzuki- Trotter ap- 
proximation and used local updates 0, Q • It has been 
replaced by the loop algorithms, wich use non-local up- 
dates. Using of non-local loop updates allows to decrease 
autocorrelation times by orders of magnitude |2|- Later 
the loop algorithms in continuous imaginary time have 
been developed Q- The continuous-time implementa- 
tion of the loop algorithm has eliminated errors result- 
ing from the Trotter discretization and, hence, loop al- 
gorithms have become numerically exact methods. 

Unfortunately, loop algorithms are inefficient in the 
presence of external field The origin of this slow- 
down results from the method of including external field 
into the simulations. External field is taken into account 
through the global weight, which increases with field in- 
crease. To construct efficient algorithm one should takes 
into account external field locally, in the loop construc- 
tion. For the first time this idea was implemented in the 
framework of the worm algorithm [f| ■ 

Both worm and loop algorithms work directly in con- 
tinuous imaginary time. At the same time there is a 
numerically exact quantum Monte Carlo method, which 
works in the discrete basis. It is Stochastic Series Expan- 
sion (SSE) method. SSE algorithm is based on power 
series expansion of a partition function. Initially SSE 
method was developed with local updates Ml. Later the 
algorithm with loop updates was proposed . Applying 
loop updates for SSE method has same favorable conse- 
quence as for world line algorithms and SSE method has 
become powerful tool for exploring quantum many-body 
systems. Recently Sandvik and Syljuasen introduce the 
concept of directed loops in stochastic series expansion, 
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which allows to perform simulation in wide range of ex- 
ternal fields @- 

Last years loop algorithms and SSE algorithm have 
been used for exploring of different quantum systems. 
Investigations of quantum spins [l£J, [!__, HtL FlI , bosons 
and one-dimensional fermion systems |12| have been 
performed. However, at the present moment investiga- 
tions of hard-core bosons and spin 5 = 1/2 systems are 
predominant in literature. 

The authors of the papers |13( investigate spin systems 
with spin S > 1/2 by loop algorithms. But they do not 
take into account external field. And they use the spin- 
split representation, i.e. replace the original spin opera- 
tors by the sum of 25 Pauli operators. Such representa- 
tion is uncomfortably because it requires extra memory 
resources and it cannot be applied directly for soft-core 
bosons. 

P. Henelius et al. have studied ferromagnetic Heisen- 
begr model with spin up to S = 2 in the wide range of 
external field by using of the SSE algorithm 0. Our 
calculations indicate that the standard SSE algorithm is 
quite effective in the case of ferromagnetic Heisenberg 
model but for simulation Heisenberg antiferromagnet it 
is necessary to increase efficiency of the algorithm. 

Until now we do not know about simulations of soft- 
core bosons by the loop or SSE algorithms. Very re- 
cently Kawashima et al. develop method for free soft- 
core bosons based on the mapping of bosonic models to 
the spin models 16] . For simulation of spin system they 
use coarse-grained loop algorithm with the spin-split rep- 
resentation. Unfortunately, the authors do not give any 
quantitative characteristics of their algorithm efficiency. 

In the present work we propose universal algorithm 
based on SSE method, which allows to investigate both 
spin systems with arbitrary spin in the presence of ex- 
ternal field and systems of interacting soft-core bosons 
in the presence of chemical potential. Also we develop 
simple procedure, which allows to increase efficiency of 
the SSE algorithm in the general case. 



II. 



THE ALGORITHM 



During the construction of the algorithm we follow the 
ideas of the work Q, therefore we do not describe SSE 
method in details but briefly outline it. 

Let us to consider Heisenberg model in the case of ar- 
bitrary spin 5, in the presence of external longitudinal 
field h 
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FIG. 1: An example of different vertices, a) In the case of 
hard-core bosons or S — 1/2 Heisenberg model, b) In the case 
of the spin-split representation for the 5 = 3/2 Heisenberg 
model, c) Vertex b), in the filling number representation. 
"1" is identified with the spin projection 5 Z = — 1/2, "2" is 
identified with the spin projection S z = 1/2. For the Bose 
Hubbard model one can identify " 1" with one boson per site, 
"2" with two bosons per site. 



where denotes summation over the pairs of nearest- 
neighbor sites. Following to the ideas of the SSE method, 
we rewrite the Hamiltonians (lllL'll as a sum over diagonal 
and off-diagonal bond operators 



<i,j> 



(3) 



where minus corresponds to antiferromagnet, plus corre- 
sponds to ferromagnet and Hubbard model (for Hubbard 
model J corresponds to t). In the case of the Heisenberg 
model the operators are 
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and, correspondingly, in the case of the Bose Hubbard 
model the operators are 
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One should guarantee non-negativity of all matrix ele- 
ments of the operators (|4I5|I by appropriate choosing of 
constants C. 

The SSE algorithm is based on the series expansion 
of the partition function Z with respect to inverse tem- 
perature (3 powers. To simplify Monte Carlo simulation, 
Sandvik et al. H, propose to introduce unit operators 
/ and cut off the expansion at n — L power. It should 
be point out, that unit operators can be distributed in 
different ways. So we obtain the formula for the partition 
function 



Z = EE 

« {Sl} 



{jp) n (L-n)\ 



k=l 
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where 7 denotes the operator type - unit, diagonal, non- 
diagonal, Sl is a sequence of operator indices and n is 
the number of non-unit operators in Sl ■ 

The Monte Carlo simulation is carried out with diag- 
onal and loop updates. The simulation starts with an 
arbitrary state \a > and operator string Sl containing 
only unit operators. During the diagonal update one at- 
tempts to interchange diagonal and unit operators with 
the probabilities 



p(f Mdh JN/3(a(p)\H^\a(p)) 



L — 
L — n - 
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JN(3(a{p)\H\f\ a {p)) 



where \a(p)) is the system state after p operators been 
applied to it, N is a number of bonds. Note that diagonal 
update changes the expansion power n by ±1. 

In the stage of loop update interchanging of diagonal 
and non-diagonal operators is carried out with the fixed 
expansion power n. At the same time system state \a > 
can be changed. 

In the case of spin S =1/2 loop update is executed in 
the following way. Non-unit operators can be represented 
as vertices with four legs (Fig.^J). One of the n vertices 
is selected and one of its four legs is selected at random. 
After that exit leg of the vertex is selected according 
with appropriate probabilities and the spins at both the 
entrance and exit legs are flipped. Note that the exit leg 
uniquely points to the entrance leg of the next vertex. 
The loop is constructed such way until it closed. 

At S > 1/2 spin-split representation of spin operators 
is widely used (Fig. Ob)). In this case vertex contains 
4(2 S + 1) variables, which can take the value ±1. During 
the construction of the loop spins at the entrance and the 
exit legs are flipped. But now loop propagates through 
the vertices with 4(25 + 1) legs. And therefore a num- 
ber of possible loop pathes increases rapidly with spin 
increase. 
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SSE algorithm allows to refuse spin-split representa- 
tion and to apply filling number representation which is 
applicable both for Heisenberg and Bose Hubbard model. 
In order to do it we use well known expressions for the 
matrix elements of corresponding operators 



(s\S , 
(n\b + \n 
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Now vertex has only four legs at arbitrary spin or at ar- 
bitrary maximum filling number for bosons (Fig. ^.j). 
However, variables connected with legs take values 
—S,...,S for spins or 0, ■■■,n ma x for bosons. Therefore 
during the construction of the loop we cannot use only 
flip of states at entrance and exit legs. So we introduce 
increasing and decreasing processes. To avoid discon- 
tinuous loop pathes during the construction of loops we 
use a simple rule: if state at the exit leg is decreased 
(increased) then at the entrance leg of the next vertex 
decreasing (increasing) process will be chosen. 



III. OPTIMIZATION OF THE ALGORITHM 

Recently Sandvik and Syljuasen Q showed that in or- 
der to fulfil detailed balance for loop update one should 
to solve the set of equations 
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(9) 



where Wi are matrix elements of the operators I|4I5[) . 
and Q>%j Eire all allowed processes. For example an de- 
notes bounce process, which does not change matrix el- 
ement Wi, and aij denotes process which transforms Wi 
to Wj. It should be point out that all must be 
nonnegative and because of detailed balance principle 
fly = aji. From atj one can obtain probabilities of all 
processes P(Wi — > Wj) = o^/Wi and correspondingly 
P{W 3 - Wi) = aij/Wj. 

We found that in the case of arbitrary spin, set of all 
processes {fty} is decomposed into closed groups con- 
taining one, three and six non-bounce processes. The 
group with one non-bounce process is described by equa- 
tions set with two equations, and groups with three and 
six non-bounce processes are described by equations sets 
with three and four equations correspondingly. So the 
equations set (0 is decomposed into sets consisting of 
two, three and four equations. Relations between num- 
ber of various groups is different at different values of 
spin. For example in the case of 5 — 1/2 there are only 
groups containing three non-bounce processes. However 
at 5 — 1 groups containing three and six non-bounce pro- 
cesses appear. And number of such groups grows with 
increase of spin until spin value becomes S = 5/2. At 
5 = 5/2 part of groups with one non-bounce process is 
4/15, with three non-bounce processes is 3/15 and with 
six non-bounce processes is 8/15. At S > 5/2 the re- 
lations between number of groups are the same as for 
5 = 5/2. 
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FIG. 2: Upper plot - bounce probabilities vs external field in 
the S — 5/2 antiferromagnetic Heisenberg model at N s = 16 
and j3 — 10. Coupling constant J = 1.0. Dark circles corre- 
spond to the optimized algorithm, open circles correspond to 
the heat-bath algorithm. Lower plot - integrated autocorrela- 
tion times for the magnetization vs external field in the cases 
of optimized and heat-bath algorithms. 



It is obvious that there is particular non-negative so- 
lution of the equations set I©- It is so-called heat-bath 
solution. 



WjWj 
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In the denominator sum is over all matrix elements be- 
longing to the group. Unfortunately heat-bath solution 
gives rise to the inefficient algorithm since all bounce pro- 
cesses an are nonzero. In order to increase efficiency of 
algorithm, one should to exclude bounce processes. Let 
us to do it for different types of groups. 

In the case of the group with one non-bounce process 
corresponding set of equations is 



W\ = an + oi2 
W 2 = a 2 2 + O21. 



(11) 



So we can always exclude one of bounce processes by 
choosing 012 = Wbjfflii = W\ — ^2,022 = if W\ > W2 
and 012 = Wi,a22 = W2 — W\,a\x = otherwise. It is 
obvious that if W\ = W% bounce processes are absent at 
all. 

Sandvik and Syljuasen for 5 = 1/2 Heisenberg model 
have analysed analytically groups with three non-bounce 
processes 0, which are described by the equations set 



W\ = an + oi2 + Oi3 
W2 = a 2 2 + 021 + a 2 3 
W 3 = a 33 + a 3 i + a 32 . 



(12) 
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FIG. 3: Integrated autocorrelation times for the magnetiza- 
tion and energy vs external field in ferromagnetic (upper plot) 
and antiferromagnetic (lower plot) Heisenberg model with dif- 
ferent spin S at N s = 16 and (5 — 10. Coupling constant is 
J = 1.0. 



FIG. 4: Integrated autocorrelation times for the mean number 
of bosons and energy vs chemical potential in the Bose Hub- 
bard model with different maximum site filling at N s — 16 
and (3 = 10. Hopping constant is t = 1.0, U = 0.5, V = 0.5. 



They proposed different solutions of the equations 
set <|12[) for various parameters of the model. It should 
be point out, that some solutions contain two bounce 
processes. At the same time for the case of arbitrary 
spin one cannot analytically analyse all allowed processes 
and obtain corresponding probabilities because number 
of processes grows rapidly as spin increase. 

We considered the equations set <|12f) in general and 
concluded that only one bounce process is needed at any 
Wi. And there is no need to solve equations set IfT^i 
analytically, but it is possible to use simple procedure 
for obtaining non-negative solution of Eq. I|12|) • 

First we demand all bounce processes an to be zero. 
Then solution of Eq. (|12|l takes the form 
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(13) 



(We take into account that o^ = a^.) If one of 
is negative then two others are certainly positive. So we 
need only one bounce process. Let a\ 2 < to be negative, 
then we should introduce bounce 033 in a such way that 
0,12 becomes positive and 0,13, 023 do not change the sign. 



Let W\ > W 2 , then by choosing 033 
we get new solution of Eq. 112|) 



W3-W1- W 2 /S 



ai 2 
a-13 
a 23 
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It is obvious that at any 5 > 1 solution Ijl4(l is positive. If 
W 2 > Wi one should to interchange W\ by W 2 in (|r4"|) . It 
should be point out that at S = 1 solution lfh4*|) coincides 
with some solutions proposed in Ref.0- We do not assert 
that our solution is most effective, but given procedure 
is universal and it can be applied at arbitrary spin. 

The groups with six non-bounce processes are de- 
scribed by the equations set 



Wi = an + a 12 + a 13 + a 1A 

W 2 — a 22 + a 21 + a 23 + a 2i 

W 3 = a 33 4- 031 4- 032 4- 034 

Wi — 044 4- 041 4- 042 4- 043- 



(15) 



As well in the case of group with three non-bounce 
processes, we demand an — and take into account 
dij = Ojj. Then we obtain the equations set with four 
equations and six variables, i.e. we have two free param- 
eters. Let us assume 023 = 034 = 013, then we obtain 
solution of Eq. ljT5|l 
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024 

We can guarantee positivity of terms like {W\ + W 2 — 
Wi)/2, by using procedure which we apply for the equa- 
tions set with three equations. Thus we introduce one 
bounce process. After that we obtain expressions like 
a — W 3 /6 with positive a. If latter expression is negative 
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one can add process 0,33 = W3 (1 — 1/82). And we can pro- 
vide positivity of solution (|16|) by choosing 82 sufficiently 
large. 



IV. TEST CALCULATIONS 

SSE algorithm is universal in any dimension. With in- 
crease of dimension extra bonds arise, but ideas of loop 
construction remain the same. Therefore we test the pro- 
posed scheme on ID systems. 

We calculate magnetization M for Heisenberg model, a 
mean number of bosons Nb for Bose Hubbard model, and 
energy for both models. We use well-known estimators 
9] 



1 " 

M = w^ {S!) (17) 

s i—1 

1 Ns 

s i=i 

where n is a number of non-unit operators in operator 
string and N s is a number of sites. We have checked 
our results with exact diagonalization and have found 
that relative deviation our results from exact is less then 
1(T 3 - 1CT 4 . 

It is well-known that integrated autocorrelation times 
is a quantitative measure of efficiency of a Monte Carlo 
sampling. We calculate autocorrelation times using bin- 
ing method, which described in Ref© 

First of all it is interesting to analyse influence of 
bounce processes on efficiency of the algorithm. To this 
end we calculate for Heisenberg antiferromagnet inte- 
grated autocorrelation times for magnetization by using 
heat-bath solution and optimized algorithm described in 
the previous section. We consider spin S = 5/2 because 
at this value all types of groups are present and relations 
between number of groups do not change with further 
spin increase. As shown in Fig. @, in the case of the 
optimized algorithm bounce probabilities are less then 
in the case of heat-bath solution. Accordingly autocor- 
relation times are less for the optimized algorithm. For 
other calculations reported here the optimized algorithm 
has been used. 

Fig. |J3J shows autocorrelation times for magnetization 
versus external field for ferromagnetic (upper plot) and 
antiferromagnetic (lower plot) Heisenberg model with 
different spin S. Calculations have been done for the 
chain with N s = 16 sites at f3 = 10. 

One can see some increase of autocorrelation times 
with spin increase for the antiferromagnet chain. How- 
ever it is difficult to compare efficiency of the algorithm 
at fixed temperature and different spin. Mean number of 
non-unit operators can be roughly estimated as N s J/3S 2 . 




t/U 



FIG. 5: Lower plot - integrated autocorrelation times for 
the mean number of bosons in the Bose Hubbard model at 
N s = 16 and = 10, N max — 5. Hopping constant is t = 1.0, 
V = 0.0, jj, = U. Dark circles correspond to the optimized al- 
gorithm, open circles correspond to the heat-bath algorithm, 
squares correspond to chain with iV s = 50. Middle plot - 
bounce probabilities for optimized and heat-bath algorithms. 
Upper plot - mean length of loops in units of (n) for optimized 
and heat-bath algorithms accordingly. The arrow points to 
the critical point "Mott insulator-superfluid" 



It is clear that this value grows rapidly with the spin S 
increase. We observe that the mean number of non-unit 
operators (n) ~ 100 at j3 = 10 and N s — 16 in the case 
of spin S = 1/2 whereas in the case of spin S = 5/2 
at the same conditions (n) ~ 2000. Thus simulation of 
S = 5/2 Heisenberg antiferromagnet at j3 = 10 is as hard 
as simulation of S = 1/2 Heisenberg antiferromagnet at 
(3 ~ 100. Hence the origin of autocorrelation time in- 
crease is clear and, with the other hand, it is obvious 
that the algorithm is very efficient. It should be point 
out that the algorithm works efficiently in wide range of 
external fields. 

For the ferromagnet chain we obtain good autocorre- 
lation times for magnetization in wide range of external 
fields with the exception of zero field. At zero field au- 
tocorrelation times for magnetization become very large 
(we do not show corresponding points at Fig. J3J). It is 
a known sequence of degeneracy states with spins up and 
spins down. 

Also we done calculations for Bose Hubbard model 
with interaction. As seen from Fig. J3J autocorrelation 
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times for energy is order of unity. Autocorrelation times 
for mean number of bosons grow with maximum filling 
number N rnax increase. Note that we can use any maxi- 
mum filling number N max , and for large class of problems 
the value N max ~ 5... 10 is quite enough. 

Investigation of many-body quantum system behav- 
ior near the critical points is one of interesting prob- 
lem in condenced matter physics. Kawashima et al. 
have tested SSE directed loop algorithm for 3D system 
and failed to obtain estimates for the observables near 
the critical point [lq|. It is well known that ID Bose 
Hubbard model experiences superfluid-insulator transi- 
tion at (t/U) c = 0.608, V = 0, (i = U 17]. We cal- 
culate autocorrelation times near the critical point for 
N s = 16, N s = 50 chains at /3 = 10, N max = 5. As 
seen from Fig. (J5J autocorrelation times for both op- 
timized and heat-bath algorithms are quite reasonable. 
But bounce probabilities in the case of heat-bath algo- 
rithm are very large and exceed bounce probabilities in 
the case of optimized algorithm by order of magnitude. 
Large bounce probabilities give rise to enormous loops, 
which walk around system many times until closed. Con- 
struction of such big loops takes a lot of time and simu- 
lation becomes inefficient. So we can conclude that SSE 
algorithm allows to perform simulations near the criti- 
cal point (at least near superfluid-insulator transition in 
ID), however it is desirable to exclude bounce processes. 



with arbitrary spin and Bose Hubbard model with inter- 
action. With the help of filling number representation we 
create the unified code for both models. Note that from 
algorithmic point of view differences between the models 
arise only at stage of matrix elements calculation. 



We propose universal procedure for excluding bounce 
processes. It has been obtained that for groups with one 
and three processes only one bounce is needed and in the 
case of group with six processes maximum two bounces 
are needed. We found that relations between number of 
various groups arc different up to spin 5 = 5/2 (maxi- 
mum filling number N max = 5). After spin S — 5/2 the 
relations do not change. 



Calculations of integrated autocorrelation times 
demonstrate increase efficiency of the algorithm under 
bounce processes excluding. We show that the proposed 
algorithm works in wide range of external fields both for 
Heisenberg model with arbitrary spin S and for Bose 
Hubbard model with interaction. Also we found that the 
algorithm is efficient near the superfluid-insulator tran- 
sition. 



V. SUMMARY 

In conclusion it should be emphasized that the algo- 
rithm introduced here allows to explore Heisenberg model 
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